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ABSTRACT 

We compute the weak coupling expansion for the energy of the three dimen- 
sional Ising model through 48 excited bonds. We also compute the magnetization 
through 40 excited bonds. This was achieved via a recursive enumeration of states 
of fixed energy on a set of finite lattices. We use a linear combination of lattices 
with a generalization of helical boundary conditions to eliminate finite volume 
effects. 



Expansions about either infinite or vanishing couphng are a major technique for the 
study of critical properties of statistical systems and field theories. These series usually 
involve a diagrammatic analysis which becomes rapidly more complex as the order in- 
creases. Thus it would be interesting to have an automated technique for the generation 
of the relevant terms. 

Here we consider generating the low temperature or weak coupling expansion for dis- 
crete systems. Our approach does not involve explicit graphs, but relies on a recursive 
computer enumeration of configurations. We illustrate the approach on the three dimen- 
sional Ising model. 

The method uses a procedure of Binder [1] for the explicit solution of discrete models 
on small lattices. In Ref. [2] these ideas were further developed. Ref. [3] explored extracting 
the low temperature series. This paper adds further tricks to obtain additional terms. This 
extends the low temperature series to an order comparable to existing high temperature 
expansions [4]. 

We consider the Ising model on a three dimensional simple cubic lattice. On each site 
i is a spin ai taking the values ±1. The system energy is 

E=J2{l-a,a,)-Hj2a, (1) 

where the first sum is over all nearest neighbor pairs of spins, each pair being counted 
once. Temporarily we set the applied field H to zero. The partition function is the sum of 
the Boltzmann weight over all configurations 

Z = J2e-^^ (2) 

Sorting configurations by energy, we rewrite this as a sum over E. We define P{E) to be 
the number of states with a given energy E. Thus, we have 
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Z=J2 P{E)u^'^ (3) 

E=0 

where A^ is the number of sites and u = e~^^. 



We compute the coefficients P{E) exactly on small systems. We recursively assemble 
the system one site at a time. The method enables us to build up a lattice with arbitrary 
length in one of the three dimensions. Intermediate stages require an explicit enumeration 
of all exposed two dimensional slices. This effectively reduces the computational complexity 
to that of a system of one less dimension. 

The starting point is a list of all states and corresponding energies for a single trans- 
verse layer of the lattice. All spins outside this layer are frozen to the same value; that is, 
the boundary conditions in the longitudinal direction are cold. Spins are then sequentially 
freed to build up the lattice in this third direction. We store the exact number of states of 
any given energy and specified exposed top layer. Storing the top layer in the bits of an 
integer /, we define p{E, I) to be this count. When a new spin or set of spins is added, we 
obtain the new counts p'{E, I) as a sum over the old counts 

p'{EJ) = J2piE-A{IJ')J'). (4) 

/' 

Here /' can differ from / only in the bits representing the newly covered spins, and A(/, /') 

is the change in energy from any newly changed bonds. For the present analysis we add 

the spins one at a time. Thus, the sum in the above equation is only over two terms, 

representing the two possible values for the newly covered spin. After the lattice is grown, 

a sum over the top layers gives the resulting P{E) = 'Yl,iP{E, I) We always continue the 

recursion sufficiently to avoid finite size errors in the longitudinal direction. 

As the temperature goes to zero, so does the variable u. Thus Eq. (3) is the low 

temperature expansion for Z. From it, we compute the series for the average energy per 

site,(£') = 2 (w^) log(.Z'). Comparing this expectation before and after adding the last 

spin, we obtain the average energy per new site. Expanding in powers of u gives 

{E/N) = Y,ejU^ (5) 

We are interested in the coefficients Cj in the infinite volume limit. One of our primary 
results, given in Table I, is the values of these coefficients through j = 48. 

At zero temperature (/3 = oo) the only surviving states have all spins parallel. As the 
temperature increases, groups of spins can fiip in this uniform background. Enumerating 



the possible combinations gives a diagrammatic method to obtain the low temperature 
coefficients [5]. Note that any enclosed group of flipped spins always involves an even 
number of excited bonds. Thus the expansion only contains even powers of u. We use 
a combination of periodic transverse and cold longitudinal boundary conditions to ensure 
this remains true on our flnite systems. 

On a periodic lattice of size nxnxn, the order to which the weak coupling expansion 
for (E/N) will agree with the inflnite volume limit is 4n — 2. A line of n flipped spins can 
wrap around the lattice and have energy 4n rather than the 4n + 2 it would have in inflnite 
space. This order can be increased via boundary conditions to require more spins to be 
flipped to wrap around the lattice. Ref. [3] showed a version of helical boundaries whereby 
an n by n transverse slice is mimicked with only [(n^ + l)/2] sites. Here we extend this 
idea to include the helicity into the longitudinal direction. 

We build our lattices one site at a time; so, it is natural to imagine the sites lying 
in a line. We do not, however, consider sequential sites as nearest neighbors. Instead, 
we introduce three integer parameters {hx,hy,hz) representing the distance along the line 
to the nearest neighbor in the corresponding x, y, or z direction. Labeling sites in the 
sequence by their ordinal number i, the nearest neighbors of site i are at i ± h^, i ^ hy 
and i ± hz- For convenience, assume hx < hy < hz- With this convention, all sites more 
than hz steps back in the chain are covered. Thus the recursion only requires us to keep 
explicit track of the hz "exposed" spins at the end of our chain. 

A minimal closed loop on such a lattice consists of a number of steps such that 
nxhx + nyhy + Uzhz = 0, where n^ represents the number of steps in the zth direction. The 
length of such a loop is n = jn^; | + |ny | + jn^ | . On an inflnite cubic lattice the only solution to 
these equations is the trivial case Ui = 0. On a flnite lattice, any other solution represents 
a flnite size correction. Flipping a chain of spins along such a closed path generates a 
state with 4n excited bonds, and creates a potential error in the series at that order. For 
a simple example, {hx, hy, hz) = (19, 21, 24) with (n^, Uy, Uz) = (3, —5, 2) gives a minimal 
loop of length 10 and will give the series to the same order as a 10^ lattice. 

Given parameters {hx,hy,hz), it is straightforward to enumerate the minimal closed 
paths. A different set of hi corresponds to a different set of such paths. However, the 



contribution to the coefficients e^ from a particular path is, by symmetry, independent of 
any permutations or sign changes in the numbers {rix, Uy, Uz). This allows us to combine 
results from various size lattices to cancel the contributions from particular closed loops. 
For example, consider loops of length 9. The (16,18,21) lattice has a minimal such loop 
with steps n = (3, 2, —4), the (16, 17, 21) lattice has closed loops with steps (1, 4, —4) and 
(5,-1,-3), the (13,18,20) lattice has a closed loops with (2,3,-4) and (4,-4,1), and 
finally the (14, 17, 19) system has the loops (3, 2, —4) and (5, —3, —1). If we combine the 
coefficients e^ as obtained from these lattices with weights (2, 1, —1, —1) respectively, then 
all errors from the loops of length 9 cancel out. This gives the series to the same order as a 
lattice with the smallest loop having length 10, which otherwise requires at least 24 sites. 

This procedure extends to cancel further loops. For our calculation, we assembled two 
lists of 10 lattices and the relative weights for combining them to cancel all loops of length 
less than 13. We ran both combinations as a check on the error cancellations. The first set 
involved lattices with hz up to 23, and took about a day on an IBM RS6000 workstation. 
The second set, involving up to hz = 24, used about half a day on a 32K Connection 
Machine and about a day of Cray-YMP time. 

Each lattice used had a minimal loop of length at least 9. While a loop of length 13 
has 52 excited bonds, we have a potential error at order 50 because of the possiblity of a 
more complex loop wrapping around the lattice simultaneously in a length 9 and a length 
10 direction. The minimal energy of such a possibility is 50 excited bonds. This is the 
limit on the order of the series presented here. 

During the recursive construction, each new count is a trivial sum of just two terms, 
representing the two possibilities for the covered spin. On the other hand, we must store 
counts for all energies up to the maximum order desired as well as for all possible values of 
the top hz spins of our helical lattice. Thus, the primary computational problem is stor- 
age. To substantially reduce these demands, we performed the calculations modulo small 
integers so that at intermediate stages the counts could be stored in one byte each. This 
gives the final coefficients modulo the given integers. After multiple passes using mutually 
prime values for these modulos, we use the Chinese remainder theorem to reconstruct the 
final series. 



From our results we constructed the series for the ratio 

As the first three ei vanish, this ratio is determined through order u^"^. te should have a 
zero at the critical point, with the slope at this zero equal to 2/ a where a is the specific 
heat exponent. 

The ratio test showed that the first singularity for the E series is unphysical and occurs 
near v? = —1/3. We therefore made a conformal transform to new variables defined by 
z = 3tt^/(l + 3tt^) to map the interval u'^ = [—1/3,0] to ^ = [— oo, 0] and the physical 
interval u'^ = [0, oo] to ^ = [0, 1]. We then did a Fade analysis in the variable z. The 
results of these are shown in Fig. 1 where we plot a few stable Fade series for te in the 
vicinity of the expected singularity in p. There is a clear zero near /3 = 0.22, in good 
agreement with the Monte Carlo studies [6] which give Pc = 0.22165.... The average slope 
of the various curves in Fig. 1 gives a = 0.22, which is about twice the accepted value for 
this exponent. The small value for this quantity makes its accurate determination difficult. 

Extending these results to include the magnetic term in Eq. (1), we augmented the 
counting to keep track of the number of flipped spins as well as excited bonds. This 
increases memory demands, so we reduced the highest energy to 20 excited bonds, and 
worked on a combination of smaller lattices with hz up to 19 to cancel closed loops of 
length 8 through 10. Assuming a spin up background, we write 

^(l-a) = 5^c,,V^A^- (7) 

where A = exp(— 2/3if). The coefficients Cij through 20 excited bonds are given in Table 
II. 

Summing the numbers in Table II over rows gives the expansion in -u^ for the mag- 
netization at zero applied field. In Fig. 2 we show several Fade approximants for the 
ratio 

-^ = ^ (8) 

Ou 

in the vicinity of the critical point. Before making these approximants, we made the same 
change of variables as used for Fig. (1). These give an estimate for fie = 0.222 and the 
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exponent (3 = 0.31, where (3 is defined by (a) oc (/3 — /9c)^ in the critical region. These 
numbers are in reasonable agreement with the accepted values. 

The method presented here should easily generalize to other discrete systems. The 
helical lattices used, as well as the combinations to cancel out finite size errors, are in- 
dependent of the Ising nature of the spins. It is straightforward to introduce additional 
couplings, although this will increase memory needs. Some interesting possibilities for 
futher exploration are gauge. Potts, and coupled gauge-spin models in various dimensions. 
Changing boundary conditions should enable the study of interface properties. A direct 
application of these counting methods to the high-temperature or strong-coupling limit 
may also be quite useful. In Ref. [7] similar recursive methods were suggested as a means 
to study many fermion systems. A particularly challenging problem is the extension of 
these ideas to theories with continuous spins. 
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Figure Captions 

Fig. 1. The ratio te defined in Eq. (6) in the vicinity of the Ising critical point. The 
series expansion for this quantity was Pade approximated in ^ = 3^^/(1 + 3w^) as the ratio 
of two polynomials, and the curves are labeled by the highest power of z appearing in the 
numerator. 

Fig. 2. The same as Fig. 1 but now for the ratio r^ in Eq. (8). 



Table I. The low temperature expansion coefficients for the average energy per unit volume. 

i Ci 



2 

4 

6 12 

8 

10 60 

12 -84 

14 420 

16 -1,056 

18 3,756 

20 -11,220 

22 37,356 

24 -118,164 

26 389,220 

28 -1,261,932 

30 4,163,592 

32 -13,680,288 

34 45,339,000 

36 -150,244,860 

38 500,333,916 

40 -1,668,189,060 

42 5,579,763,432 

44 -18,692,075,820 

46 62,762,602,860 

48 -211,062,133,044 



Table 11. 


Coefficients Qj 


for the 


expansion of the 


mag] 


aetization. 


Here \{l - a) = 


E^,J C^JU^ 


''A-^ where u = e 


-2/3 and A = e- 


2/3H Unhsted coefficients for i < 20 all van- 
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-1,101,660 
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4,722,564 
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-7,580,660 
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8,680,245 


51,142,152 
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-4,397,652 


-130,897,242 
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216 
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1,143 
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49,896 




3,960 




264 




36 


17 


360,450 




41,310 




7,26C 


) 





18 


547,236 




672,670 




73,216 


12,960 
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-12,320,586 


2,368,080 


773,025 


138,744 


20 


-35,804,700 


-6,147,840 


6,632,208 


1,220,220 
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